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| Abstract 

An iterative scheme for the Dynamical Systems Method (DSM) is 
given such that one does not have to solve the Cauchy problem occuring 
in the application of the DSM for solving ill-conditioned linear algebraic 
systems. The novelty of the algorithm is that the algorithm does not 
■ have to find the regularization parameter a by solving a nonlinear 

equation. Numerical experiments show that DSM competes favorably 
1****" ■ with the Variational Regularization. 
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^ ; 1 Introduction 

> ■ 

• i-H , 

^ , Consider a linear operator equation of the form 

F(u) = Au-f = 0, u£ H, (1) 

where if is a Hilbert space and A is a linear operator in H which is not 
necessarily bounded but closed and densely defined. To solve this equation 
we apply the Dynamical Systems Method (DSM) introduced in [6]: 

u' = -u+ (T + a(s))- l A*f, u(0)=u , (2) 

where T := A* A and a(t) > is a nonincreasing function such that a(t) — > 
as t — > oo. The unique solution to ([2]) is given by 



u{t) = u e- t + e- t [ e s {T + a{s))- 1 A*fds. (:!) 
J o 
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The DSM consists of solving problem ([2]) with a chosen a(t) and uq 
and finding a stopping time tg so that u(tg) approximates the solution y 
to problem (TjQ) of minimal norm. Different choices of a(t) generate different 
methods of solving equation ([1]) . These methods have different accuracy and 
different computation time. Thus, in order to get an efficient implementation 
of the DSM, we need to study the choice of a(t) and of the stopping time 
tg. Since the solution to (H|) can be presented in the form of an integral, 
the question arises: how can one compute the integral efficiently? The 
integrand of the solution is used also in the Variational Regularization (VR) 
method. The choice of the stopping time tg will be done by a discrepancy- 
type principle for DSM. However, choosing a(t) so that the method will be 
accurate and the computation time is small is not a trivial task. 

This paper deals with the following questions: 

1. How can one choose a(t) so that the DSM is fast and accurate? 

2. Does the DSM compete favorably with the VR in computation time? 

3. Is the DSM comparable with the VR in accuracy? 

2 Construction of method 

2.1 An iterative scheme 

Let us discuss a choice of a(t) which allows one to solve problem (|2|) or to 
calculate the integral ([3]) without using any numerical method for solving 
initial- value problem for ordinary differential equations (ODE). In fact, using 
a monotonically decreasing a(t) with one of the best numerical methods for 
nonstiff ODE, such as DOPRI45, is more expensive computationally than 
using a step function a(t), approximating a(t), but brings no improvement 
in the accuracy of the solution to our problems compared to the numerical 
solution of our problems given in Section [3.1.21 

Necessary conditions for the function a(t) are: a(s) is a nonincreasing 
function and lim s _ >00 a(s) = (see [6]). Thus, our choice of a(t) must satisfy 
these conditions. Consider a step function a(t), approximating a(t), defined 
as follows: 

a{t) = a n , t n <t < t n+ i, 
the number t n are chosen later. For this a(t), u n = u(t n ) can be computed 
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by the formula: 



n 

i=l 

This leads to the following iterative formula: 

w n+ i = e- h "-u n + (l-e~ /l ")(r + a n ) _1 ^*/ 5 , K = t n+1 -t n . (4) 

Thus, ti ra can be obtained iteratively if uq , a(i) and t n are known. 
The questions are: 

1. For a given a(t), how can we choose t n or /i n so that the DSM works 
efficiently? 

2. With a n = a(t n ) where a(t) is a continuous function, does the iterative 
scheme compete favorably with the DSM version in which u(t) is solved 
by some numerical methods such as Runge-Kutta methods using a(i)? 

In our experiments, a n = a{t n ) where a(t) = jfi where ao > is a 
constant which will be chosen later, as suggested in [6], with t n chosen so 
that t n+ \ — t n = h n , h n = q n , where 1 < q < 2. For this choice, if q > 1 
then the solution u n at the n-th step depends mainly on (T + a n ) 1 A* fg 
since e~ hn is very small when n is large. 

Note that a n decays exponentially fast when n — * oo if q > 1. A question 
arises: how does one choose q so that the method is fast and accurate? This 
question will be discussed in Section 3. 

ALGORITHM Q] demonstrates the use of the iterative formula (j4|) and a 
relaxed discrepancy principle described below for finding u given ao, A, fs 
and 5. 

In order to improve the speed of the algorithm, we use a relaxed discrep- 
ancy principle: at each iteration one checks if 

0.95 < \\Aun-fs || < 1.0015. (5) 

As we shall see later, ao is chosen so that the condition ([7|) (see below) is 
satisfied. Thus, if no = T~ ) 1 A*fs, where T a := T + a, then 5 < \\Auo — fs\\. 
Let t n be the first time such that \\Au n — fg\\ < 1.0015. If © is satisfied, 
then one stops calculations. If \\Au n — f$\\ < 0.95, then one takes a smaller 
step-size and recomputes u n . If this happens, we do not increase h n , that 
is, we do not multiply h n by q in the following steps. One repeats this 
procedure until condition ([6]) is satisfied. 
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Algorithm 1 : BSM(A,f s ,S) 
q:=2; 

g 5 :=A*f 5 ; T:=A*A; 

itermax := 30 ; u = (T + ao)^ 1 gs ; 

i:=0; t = l; h:=l; halve := ; 

while (1.0015 < \\Au — fg\\) and {i < itermax) do 

i:=i + l; t = t + h; a = ao/t; 

v := (T + al^gs; 

u = e~ h u + (1 — e~ h )v ; 

if 0.95 < \\Au - f s \\ then 
u := u ; 

if halve = then h := hq ; end ; 
elseif 

t := t — h; h := /i/2 ; halve = 1 ; 
endif 
endwhile 



In order to improve the speed of the algorithm, we use a relaxed discrep- 
ancy principle: at each iteration one checks if 

0.95 < \\Au n - f 5 \\ < 1.0015. (6) 

As we shall see later, ao is chosen so that the condition © (see below) is 
satisfied. Thus, if uq = T~ ) 1 A*fs, where T a := T + a, then 5 < \\Auo — fs\\. 
Let t n be the first time such that \\Au n — fs\\ < 1.0015. If (jGj) is satisfied, 
then one stops calculations. If \\Au n — f$\\ < 0.95, then one takes a smaller 
step-size and recomputes u n . If this happens, we do not increase h n , that 
is, we do not multiply h n by q in the following steps. One repeats this 
procedure until condition ([6]) is satisfied. 

2.2 On the choice of ao 

From numerical experiments with ill-conditioned linear algebraic systems 
(las) of the form Au = fs, it follows that the regularization parameter om, 
obtained from the discrepancy principle ||An ajlf — fs\\ = 5, where u au = 
T~ A 1 I A*fs, is often close to the optimal value a op , i.e., the value minimizing 
the quantity: 

IK™ - y\\ = inf |K - y\\, u a = T^Afs. 

y a 
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The letter M in clm stands for Morozov, who suggested to choose c = 1 in 
the disrepancy principle. 

If ao is chosen smaller than a op , the method may converge poorly. Since 
aM is close to a op , only those a for which \\AT~ 1 A* fs — fs\\ = cS with c 
'close' to 1 yield accurate approximations to the solution y. Also, if ao 
is chosen much greater than a op , then the information obtained from the 
starting steps of the iterative process (HJ) is not valuable because when ao 
is far from a op , the error \\u ao — y\\ is much bigger than \\u aop — y\\. If ao is 
much bigger than a op , a lot of time will be spent until a(t n ) becomes close 
to a op . In order to increase the speed of computation, ao should be chosen 
so that it is close to a op and greater than a op . Since a op is not known and is 
often close to aM, we choose ao from the condition: 

5 < \\Au ao - f 5 \\ < 26. (7) 

For this choice, ao is 'close' to and greater than aM- Since there are many 
ao satisfying this condition, it is not difficult to find one of them. 

In the implementation of the VR using discrepancy principle with Mo- 
rozov's suggestion c = 1, if one wants to use the Newton method for finding 
the regularization parameter, one also has to choose the starting value ao so 
that the iteration process converges, because the Newton method, in gen- 
eral, converges only locally. If this value is close to and greater than aM, it 
can also be used as the initial value of ao = a(t)\ t= o for the DSM. 

In our numerical experiments, with a guess ao = | max \i(A*A)5 re i for 
a(0), we find ao such that 5 < \\Au ao — fs\\ < 25. Here, 5 re i stands for the 
relative error, i.e., 5 re i = JLr. The factor | is introduced here in order to 
reduce the cost for finding ao, because ao, which satisfies ([7]), is often less 
than max \((A* A)5 re i. The idea for this choice is based on the fact that the 
spectrum of the matrix max xl(A*A) ^*^ * s contained in [0, 1]. 

Note that ones has 

oP|| 2 

aM ~ wTt- 

Indeed, 

\\fs\\-8=\\h\\-\\Au aM -fs\\<\\Au au \\. 

Since A*Au aM + a M u aM = A*f s , one has a M Au aM = AA*(f s - Au au ). 
Thus, 

\\fs\\-5\\ < \\Au aM \\ = J-\\AA*(f 5 -Au aM )\\ < ^o. 

aM aM 
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Similar estimate one can find in (5J p. 53], where ao = ||^||_^ * s suggested as 
a starting value for Newton's method to determine ajf on the basis that it 
is an upper bound for a,M- Note that [fM~5 ~ ^rezll^ll 2 = ToaaxXi(A*A)6 re i. 
However, in practice Newton's method does not necessarily converge with 
this starting value. If this happens, a smaller starting value ai := is used 
to restart the Newton's method. 

In general, our initial choice for ao may not satisfy ©. Iterations for 
finding ao to satisfy (JT]) are done as follows: 

1. If ^ = c > 3, then one takes a\ := 2 (c-i) as ^ ne nex ^ S uess 

and checks if the condition (|7|) is satisfied. If 2 < c < 3 then one takes 
ai := a /3. 

2. If ^ Aua( > s ^ = c < 1, then a\ := 3ao is used as the next guess. 

3. After ao is updated, one checks if ([7]) is satisfied. If ([7]) is not satisfied, 
one repeats steps 1 and 2 until one finds ao satisfying condition ([7]) 
(see ALGORITHM©. 



Algorithm 2 : fmd-a 

111 A 1 1 2 r 
°0 : = 3 H^MI re i 

c:= \\Au ao - fs\\/5; 
while (2 < c) or (c < 1) do 
if 3 < c then 

a := 0.5a /(c - 1) ; 
elseif (2 < c < 3) then 

a := a /3 ; 
else 

ao := 3ao ; 
end 

u ao := (A*A + a )- 1 A*f s ; 
c := \\Au ao - fsW/8 ■ 
endwhile 



The above strategy is based on the fact that the function 
<P(a) = \\A(T + a)- 1 A*fs-fs\\ 
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is a monotonically decreasing function of a, a > 0. In looking for ao, sat- 
isfying ((TJ), when our guess ao 3> aju > or ||^4u ao — /$|| S> 5, one uses an 
approximation 

At \ ~ At \ , c >( a o) - 

« 0(a o ) + [x - ao) 

ao — ajw 

« 0(a o ) + (x - a ) =: ^ J . 

a 

Note that </>(ao) and ao are known. We are looking for x such that 5 < 
<p(x) < 25. Thus, if a\ is such that 5 < <p(ai) < 26 and if 25 < 4>(ao), then 

((p(a ) - 25) Q ° — - < a - ai < (<p(a ) - 5)- 



A (a ) -5 <p(a ) - 5' 

Hence, we choose a± such that 

a 



ao - a\ = (0(oo) - 1-55)- 
0.5o 



0(a o ) - 5' 

so 

ai = a 



</>(a ) - 5 

Although this is a very rough approximation, it works well in practice. It 
often takes 1 to 3 steps to get an ao satisfying ([7]). That is why we have 
a factor in the first case. Overall, it is easier to look for ao satisfying 
([7|) than to look for ao for which the Newton's method converges. Indeed, 
the Newton's scheme for solving om does not necessarily converge with ao 
found from condition ([7]). 



3 Numerical experiments 

In this section, we compare DSM with VRj and VR n . In all methods, we 
begin with the guess ao = |||^4|| 2 (5 re / and use the ALGORITHM [2] to find ao 
satisfying condition ([TJ. In our experiments, the computation cost for this 
step is very low. Indeed, it only takes 1 or 2 iterations to get ao- By VRj we 
denote the VR obtained by using a = ao, the intial value for a{t) in DSM, 
and by VR n we denote the VR with a = clm , found from the VR discrepancy 
principle with c = 1 by using Quasi-Newton's method with the initial guess 
a = a,Q. Quasi-Newton's method is chosen instead of Newton's method in 
order to reduce the computation cost. In all experiments we compare these 
methods in accuracy and with respect to the parameter Nu nso i, which is the 
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number of times for solving the linear system T a u = A*fs for u. Note that 
solving these linear systems is the main cost in these methods. 

In this section, besides comparing the DSM with the VR for linear al- 
gebraic systems with Hilbert matrices, we also carry out experiments with 
other linear algebraic systems given in the Regularization package in [3J. 
These linear systems are obtained as a part of numerical solutions to some 
integral equations. Here, we only focus on the numerical methods for solving 
linear algebraic systems, not on solving these integral equations. Therefore, 
we use these linear algebraic systems to test our methods for solving stably 
these systems. 

3.1 Linear algebraic systems with Hilbert matrices 

Consider a linear algebraic system 



H n u = fs, 



(8) 



where 



fs = f + e, f = H n x, H n 



l 

n+l 



1 

n+l 



1 

2n-l. 



E WL n is a random normally distributed vector such that ||e||2 < 
2- The Hilbert matrix H n is well-known for having a very large 



and e 

SrelWf 

condition number when n is large. If n is sufficiently large, the system is 
severely ill-conditioned. 



3.1.1 The condition numbers of Hilbert matrices 

It is impossible to calculate the condition number of H n by computing the 
ratio of the largest and the smallest eigenvalues of H n because for large n 
the smallest eigenvalue of H n is smaller than 10 -16 . Note that singular val- 
ues of H n are its eigenvalues since H n is self adjoint and positive definite. 
Due to the limitation of machine precision, every value smaller than 10~ 16 is 
understood as 0. That is why if we use the function cond provided by MAT- 
LAB, the condition number of H n for n > 20 is about 10 16 x max \ \i(H n )\. 
Since the largest eigenvalue of H n grows very slowly, the condition numbers 
of H n for n > 20 are all about 10 20 , while, in fact, the condition number of 
H W q computed by the formula, given below, is about 10 150 (see Table QJ. In 
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general, computing condition numbers of strongly ill-conditioned matrices 
is an open problem. The function cond, provided by MATLAB, is not al- 
ways reliable for computing the condition number of ill-condition matrices. 
Fortunately, there is an analytic formula for the inverse of H n . Indeed, one 
has (see [2]) H~ l = (hij)™ j=1 , where 

Thus, the condition number of the Hilbert matrix can be computed by the 
formula: 

cond(H n ) = WHnWWH^W. 

Here cond(H n ) stands for the condition number of the Hilbert matrix H n and 
\\H n \\ and \\H~ || are the largest eigenvalues of H n and H~ l , respectively. 
Although MATLAB can not compute values less than 10 -16 , it can compute 
values up to 10 200 . Therefore, it can compute H-f^ 1 )! for n up to 120. 
In MATLAB, the matrices H n and H^ 1 can be obtained by the syntax: 
H n = hilb(n) and H~ l = invhilb(n), respectively. 

The condition numbers of Hilbert matrices, computed by the above for- 
mula, are given in Table [TJ 

Table 1: The condition number of Hilbert matrices. 



n 


20 


40 


60 


80 


100 120 


cond(H n ) 


2.5 x 10^ 


7.7 x 10 bs 


2.7 x 10 sy 


9.9 x 10 liy 3i 


5 x 10 lbU 1.5 x 10 i8i 



From Table [U one can see that the computed condition numbers of the 
Hilbert matrix grow very fast as n grows. 

3.1.2 Continuous a(t) compared to the step function d(t) 

In this section, we compare the DSM, which is implemented by solving the 
Cauchy problem ([2]) with a(t), and the iterative DSM implemented with d{t) 
approximating a(t) as described in Section [2.11 Both of them use the same 
ao which is found by ALGORITHM [21 The DSM using a numerical method 
to solve the Cauchy problem is implemented as follows: 

1. One uses the DOPRI45 method which is an embedded pair consisting 
of a Runge-Kutta (RK) method of order 5 and another RK method of 
order 4 which is used to estimate the error in order to control the step 
sizes. The DOPRI45 is an explicit method which requires 6 right-hand 
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side function evaluations at each step. Details about the coefficients 
and variable step size strategy can be found in [H [3] . Using a variable 
step size helps to choose the best step sizes and improves the speed. 

2. In solving ([2]), at the end of each step, one always checks the stopping 
rule, based on the discrepancy principle 

0.9 < \\Au s (t) - f s \\ < 1.001(5. 

If this condition is satisfied, one stops and takes the solution at the 
final step u(t n ) as the solution to the linear algebraic system. 



Table 2: Numerical results for Hilbert matrices for 5 re i = 0.01, n = 100. 







DSM 


DSM(<j = 1) 


DSM-DOPRI45 


n 


Nuns 




-Wlinsol 


ll«a-y||2 


-^Vlinsol 


l|w«— v||a 


ol 1Mb 


1Mb 


llj/lb 


10 


5 


0.1222 


10 


0.1195 


205 


0.1223 


20 


5 


0.1373 


7 


0.1537 


145 


0.1584 


30 


7 


0.0945 


20 


0.1180 


313 


0.1197 


40 


5 


0.2174 


7 


0.2278 


151 


0.2290 


50 


6 


0.1620 


14 


0.1609 


247 


0.1609 


60 


6 


0.1456 


16 


0.1478 


253 


0.1480 


70 


6 


0.1436 


13 


0.1543 


229 


0.1554 


80 


6 


0.1778 


10 


0.1969 


181 


0.1963 


90 


6 


0.1531 


13 


0.1535 


307 


0.1547 


100 


7 


0.1400 


23 


0.1522 


355 


0.1481 



The DSM version implemented with the DOPRI45 method is denoted 
DSM-DOPRI45 while the other iterative version of DSM is denoted just by 
DSM. 

Table [2] presents the numerical results with Hilbert matrices H n ob- 
tained by two versions of the DSM for n = 10, 20, 100, 5 re i = 0.01, 

x = (xi, x n ) T , Xi = J 2^-71". From Tabled as well as other numerical ex- 
periments, we found out that the accuracy obtained by the DSM-DOPRI45 
is worse than that of the iterative DSM. Moreover, the computation time 
for the DSM-DOPRI45 is much greater than that for the iterative DSM. 
Also, using h =const or q = 1 does not give more accurate solutions while 
requires more computation time. 

The conclusion from this experiment as well as from other experiments 
is that the DSM with q = 2 is much faster and often gives better results 
than the DSM with q = 1 and the DSM-DOPRI45. Therefore, we choose 
the iterative DSM with q = 2 to compare with the VR n method. 
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3.1.3 DSM compared to VR 



In this section, we test three methods: the DSM, the VRj and the VR n on 
linear algebraic systems with Hilbert matrices. The first linear system is 
obtained by taking -ffioo and x = (xi, ...,xioo) T , where Xi = (j^) 2 - For the 
second problem we just change X{ to Xi = sin(2^y7r). Numerical results for 
these systems are shown in Figure [TJ 




20 40 60 80 100 20 40 60 80 100 



Figure 1: Plots of solutions obtained by the DSM and VR with the exact 
solution x, x = (£Cj)*™ when X{ = (2^-7r) 2 (left) and Xi = sin(2^7r) (right) 
with 5 re i = 0.02. 

Looking at Figure [H one can see that with the same guess ao, both the 
VR n and DSM give better results than those of VRj. As it can be seen from 
Figure [IJ the numerical solutions obtained by the DSM in these tests are 
slightly more accurate than those of the VR n . 

Table [3] presents results with Hilbert matrices H n for n = 10,20, 100, 

5 re i = 0.01, x = (x\, ...,x n ) T , X{ = yj 2^jjj7r. Looking at this Table it is clear 
that the results obtained by the DSM are slightly more accurate than those 
by the VR n even in the cases when the VR n requires much more work than 
the DSM. In this example, we can conclude that the DSM is better than the 
VR n in both accuracy and time of computation. 
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Table 3: Numerical results for Hilbert matrix H n for d re i = 0.01, 
10,20, ...,100. 



n 



10 
20 
30 
40 
50 
GO 
70 
80 
90 
100 



linsol 



DSM 

us — |<j[g 

vh 



0.2368 
0.1638 
0.1694 
0.1984 
0.1566 
0.1890 
0.1449 
0.1217 
0.1259 
0.1865 



VR,; 



\Ws-y\\2 

II y II 2 



0.3294 
0.3194 
0.3372 
0.3398 
0.3345 
0.3425 
0.3393 
0.3480 
0.3483 
0.2856 



linsol 



vr„ 



\U\V2 



7 
7 
11 
8 
7 
8 
11 
8 
11 
9 



0.2534 
0.1765 
0.1699 
0.2074 
0.1865 
0.1980 
0.1450 
0.1501 
0.1355 
0.1937 



3.2 A linear algebraic system related to an inverse problem 
for the heat equation 

In this section, we apply the DSM and the VR to solve a linear algebraic 
system used in the test problem heat from Regularization tools in [JJ. This 
linear algebraic system is a part of numerical solutions to an inverse prob- 
lem for the heat equation. This problem is reduced to a Volterra integral 
equation of the first kind with [0, 1] as the integration interval. The kernel 
is K(s, t) = k(s — t) with 

r 3/2 l 

^ ) = 2^ eXp( -4^ } - 

Here, we use the default value k = 1. In this test in [JJ the integral equa- 
tion is discretized by means of simple collocation and the midpoint rule 
with n points. The unique exact solution u n is constructed, and then the 
right-hand side b n is produced as b n = A n u n (see [3J). In our test, we use 
n = 10, 20, 100 and b n ^ = b n + e n , where e n is a vector containing ran- 
dom entries, normally distributed with mean 0, variance 1, and scaled so 
that ||e n || = <5 re z||& n ||. This linear system is ill-posed: the condition num- 
ber of ^4ioo obtained by using the function cond provided in MATLAB is 
1.3717 x 10 37 . As we have discussed earlier, this condition number may 
be not accurate because of the limitations of the program cond provided 
in MATLAB. However, this number shows that the corresponding linear 
algebraic system is ill-conditioned. 

Looking at the Table [Hone can see that in some situations the VR n is not 
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Table 4: Numerical results for inverse heat equation with 5 re i = 0.05, n 

m,i = mo. 



10 

20 
30 
40 
50 
GO 
70 
80 
90 
100 



linsol 



DSM 

vh 



0.2051 
0.2198 
0.3691 
0.2946 
0.2869 
0.2702 
0.2955 
0.2605 
0.2616 
0.2588 



linsol 



VRi 

\\us-yh 
vh 



0.2566 
0.4293 
0.4921 
0.4694 
0.4780 
0.4903 
0.4981 
0.4743 
0.4802 
0.4959 



VR n 

Ml II ay II 2 



6 
8 
6 
8 
7 
9 
6 

10 



0.2066 
0.2228 
0.3734 
0.2983 
0.3011 
0.2807 
0.3020 
0.2513 
0.2692 
0.2757 



as accurate as the DSM even when it takes more iterations than the DSM. 
Overall, the results obtained by the DSM are often slightly more accurate 
than those by the VR n . The time of computation of the DSM is comparable 
to that of the VR n . In some situations, the results by VR n and the VRj 
are the same although it uses 3 more iterations than does the DSM. The 
conclusion from this Table is that DSM competes favorably with the VR n 
in both accuracy and time of computation. 

Figure [2] plots numerical solutions to the inverse heat equation for S re i = 
0.05 and 5 re i = 0.02 when n = 100. From the figure we can see that the 
numerical solutions obtained by the DSM are about the same those by the 
VR n . In these examples, the time of computation of the DSM is about the 
same as that of the VR n . 

The conclusion is that the DSM competes favorably with the VR n in 
this experiment. 

3.3 A linear algebraic system for the computation of the 
second derivatives 

Let us do some numerical experiments with linear algebraic systems arising 
in a numerical experiment of computing the second derivative of a noisy 
function. 

The problem is reduced to an integral equation of the first kind. A linear 
algebraic system is obtained by a discretization of the integral equation 
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whose kernel K is Green's function 

{S,t) \t{ S -l), if 8>f 

Here s, t E [0, 1] and as the right-hand side / and the corresponding solution 
u one chooses one of the following (see [1]): 

case 1, f(s) = (s 3 — s)/6, u(t) = t, 

case 2, f(s) = e s + (1 - e)s - 1, = e\ 

o f ^ _ / (4s 3 - 3s)/24, if s < \ 

cased, /(sj - | ( _ 4s 3 + 12s2 _ 9s + 1)/24; . f g ^f, 



u(t) 



t, if t < | 
1 - t, if t > | 



Using ^4 n and ii ra in [3], the right-hand side b n = A n u n is computed. 
Again, we use n = 10, 20, 100 and b nj $ = b n + e n , where e n is a vector 
containing random entries, normally distributed with mean 0, variance 1, 
and scaled so that ||e n || = cVe/H^nll- This linear algebraic system is mildly 
ill-posed: the condition number of ^4ioo is 1.2158 x 10 . 

Numerical results for the third case is presented in Table [5j In this 
case, the results obtained by the VR n are often slightly more accurate than 



14 



those of the DSM. However, the difference between accuracy as well as the 
difference between time of computation of these methods is small. Numerical 
results obtained by these two methods are much better than those of the 
VR,. 



Table 5: Results for the deriv2 problem with 5 re i = 0.01, n = 100 case 3. 



10 
20 
30 
40 
50 
GO 
70 
80 
90 
100 



DSM 



insol 



\Ws-y\\-2 

HE 



0.0500 
0.0584 
0.0690 
0.0367 
0.0564 
0.0426 
0.0499 
0.0523 
0.0446 
0.0399 



Ni 



insol 



VR 

\us—y\\2 

lklb_ 



0.0542 
0.0708 
0.0718 
0.0454 
0.0565 
0.0452 
0.0422 
0.0516 
0.0493 
0.0415 



VR„ 



Ni 



insol 



\\us-y\\2 

\\yh 



0.0444 
0.0561 
0.0661 
0.0384 
0.0562 
0.0407 
0.0372 
0.0498 
0.0456 
0.0391 



For other cases, case 1 and case 2, numerical results obtained by the 
DSM are slightly more accurate than those by the VRj. Figure [3] plots the 
numerical solutions for these cases. The computation time of the DSM in 
these cases is about the same as or less than that of the VR n . 



The deriv2 problem with n=100, 8 re| =0.02, easel 
0.12 1 1 1 1 1 1 




20 40 60 80 100 



The deriv2 problem with n=100, 8 re| =0.02, case 2 
0.35 1 1 1 1 1 1 




Figure 3: Plots of solutions obtained by DSM, VR for the deriv2 problem 
when n = 100, 8 ret = 0.02 (left) and 5 re i = 0.02 (right). 
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The conclusion in this experiment is that the DSM competes favorably 
with the VR. Indeed, the VR„ is slightly better than the DSM in case 3 but 
slightly worse than the DSM in cases 1 and 2. 

4 Concluding remarks 

The conclusions from the above experiments are: 

1. The DSM always converges for a(t) = given that ao > a op . How- 
ever, if ao is not well chosen, then the convergence speed may be slow. 
The parameter ao should be chosen so that it is greater than and close 
to the optimal a op , i.e., the value minimizing the quantity: 

IK op - y\\ = inf IK - y\\, u a = T~ l Af 5 . 

a 

However, since a op is not known and om is often close to a op , we choose 
ao so that 

5<\\AT- 1 A*f s -f s \\<28. 

2. The DSM is sometimes faster than the VR. In general, the DSM is 
comparable with the VR n with respect to computation time. 

3. The DSM is often slightly more accurate than the VR, especially when 
5 is large. Starting with ao such that 5 < WAT^ A* f$ — fs\\ < 25, the 
DSM often requires 4 to 7 iterations, and main cost in each iteration 
consists of solving the linear system T a u = A* f$. The cost of these 
iterations is often about the same as the cost of using Newton's method 
to solve au in the VR ra . 

4. For any initial a such that 5 < \\AT~ 1 A*f s - f s \\ < 26, the DSM 
always converges to a solution which is often more accurate than that 
of the VR n . However, with the same initial ao, the VR n does not 
necessarily converge. In this case, we restart the Newton scheme to 
solve for the regularization parameter with initial guess a\ = ^ instead 
of ao- 
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